/***************************************************************************
 *   Copyright (C) 2007 by Alexandre Masselot,,,   *
 *   alex@alex-laptop   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU Library General Public License as       *
 *   published by the Free Software Foundation; either version 2 of the    *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this program; if not, write to the                 *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#include "utils/psmnistmspreader.h"
#include "msmscoring/scoringmodeldancik.h"
#include "utils/runlcmsmsreadermgf.h"
#include "utils/deltamassbuilder.h"
#include "utils/peaktheomsmsfactorybasic.h"
#include "utils/aminoacidinstancerbasic.h"
#include "utils/dbsequencefactoryfasta.h"
#include "utils/enzymetrypsin.h"
#include "utils/peptidemsmsfactory.h"

using namespace std;
using namespace bol;

#define STRINGIFY(x) XSTRINGIFY(x)
#define XSTRINGIFY(x) #x



/**
 *
 * @param argc
 * @param argv
 * @return
 */
int main(int argc, char **argv){

  //init amino acids
  AminoAcidInstancerBasic aainstancer;
  aainstancer.setMasses();
  AminoAcid::initChar2mass();

  /*
   * All the modif set and reader stuf is taken as-is from PMSNistReaderTest
  */
  string mspfile((argc>1)?argv[1]:STRINGIFY(DATADIR) "/a.msp");

  ModificationDictionary dico;
  dico.put(new Modification("Cys_CAM C 57.052"));
  dico.put(new Modification("CAM_core K 43.02504"));
  dico.put(new Modification("Oxidation_HW HW 15.999"));
  dico.put(new Modification("Oxidation_M M 15.999"));
  
  PeptideModifier peptmodifier;
  peptmodifier.setModifdico(&dico);
  peptmodifier.addModif(dico["Cys_CAM"], true);
  peptmodifier.addModif(dico["CAM_core"], true);
  peptmodifier.addModif(dico["Oxidation_HW"], true);
  peptmodifier.addModif(dico["Oxidation_M"], true);

  PSMNistMSPReader reader;
  reader.setPeptideModifier(&peptmodifier);
  reader.open(mspfile);

  //end of common block


  //tol to be used on fragments
  DeltaMassBuilder builder;
  DeltaMass *fragTol=builder.deltamassInstance("1000ppm");
  DeltaMass *precTol=builder.deltamassInstance("2Da");

  //must store a run to perform the null scanning
  RunLcmsms run;

  //fragmentation model
  PeakTheoMsmsFactoryBasic peaktheofactory;
  peaktheofactory.setSeries("y,b,b++,y++",2);

  //model declaration
  ScoringModelDancik scoringmodel;
  scoringmodel.series_add("b");
  scoringmodel.series_add("y");
  scoringmodel.init();
  scoringmodel.learn_init();

  //tempdata structure
  Spectrum<PeakTheoMsms*> theospectrum;
  SpectraAlignment<PeakTheoMsms, PeakIntensity> spalign;

  //learn the correct matches
  cout<<"pept\terr.Da\terr.ppm\n";
  while(PeptideSpectrumMatch *pmatch=reader.nextMatch()){
    peptmodifier.printPeptide(cout, pmatch->peptide());
    double d=(pmatch->peptide()->getMass()/pmatch->spectrum()->charges()[0]+1.00728-pmatch->spectrum()->moz());
    cout<<"\t"<<pmatch->peptide()->getMass()<<"\t"<<d<<"\t"<<d/pmatch->spectrum()->moz()<<"\n";
    run.msmsCollection().push_back(pmatch->spectrum());
    theospectrum.clear();
    peaktheofactory.computeTheoMasses(*pmatch->peptide(), theospectrum, pmatch->spectrum()->charges()[0], true);
    spalign.spectrumA(&theospectrum);
    spalign.spectrumB((Spectrum<PeakIntensity*>*)(pmatch->spectrum()));
    spalign.align(fragTol);
    scoringmodel.learn_pushAlternate(spalign);
  }

  //learn null matches
  run.buildIndex();
  run.msmsCollection().index()->setTolerance(precTol);
  
  //set the enzyme and a vector to store cleaved peptides
  Enzyme *enzyme=new EnzymeTrypsin();
  vector<Peptide> vcleavedpept;
    
  PeptideModifier peptmodifierNull;
  
  //the factory will be set up to produce peptide according to the given modifier (modifier cannot be modified thereafter)
  PeptideMSMSFactory factory;
  factory.setPeptideModifier(&peptmodifierNull);
  //the vector in which will be stored the msms combination
  vector<Peptide> vmsmspept;

  DBSequenceFactoryFasta fasta;
  fasta.openFile(STRINGIFY(DATADIR) "/little.fasta");

  DBSequence entry;
  int i=0;
  while(fasta.nextSequence(entry)){
    if(i++%1000==0)
      cout<<entry.ac()<<endl;
    enzyme->cleave(*entry.sequence(), vcleavedpept);
    for(vector<Peptide>::iterator it=vcleavedpept.begin(); it<vcleavedpept.end(); it++){
      vmsmspept.clear();
      factory.createPeptides(*it, vmsmspept);
      for(vector<Peptide>::iterator it2=vmsmspept.begin(); it2<vmsmspept.end(); it2++){
        vector<const SpectrumLcmsmsMassIndexAnchor*>::const_iterator first;
        vector<const SpectrumLcmsmsMassIndexAnchor*>::const_iterator last;
        run.msmsCollection().index()->getAnchors(it2->getMass(), first, last);
        for(vector<const SpectrumLcmsmsMassIndexAnchor*>::const_iterator it=first; it<last; it++){
          theospectrum.clear();
          peaktheofactory.computeTheoMasses(*it2, theospectrum, 2, true);
          spalign.spectrumA(&theospectrum);
          spalign.spectrumB((Spectrum<PeakIntensity*>*)(*it)->spectrum());
          spalign.align(fragTol);

          scoringmodel.learn_pushNull(spalign);

          //double score=scoringmodel.score(spalign);
          //cout<<"\t"<<score<<"\t"<<((*it)->spectrum())->description()<<endl;
        }
      }
    }
  }


  scoringmodel.learn_finalize();
  cout<<scoringmodel<<endl;
}
